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To preserve nonlinearity of a full order system over a parameters range of interest, we propose a simple 
modeling approach by assembling a set of piecewise local solutions, including the first-order Taylor series 
terms expanded about some sampling states. The work by Rewienski and White [1] inspired our use of 
piecewise linear local solutions. The assembly of these local approximations is accomplished by assigning 
nonlinear weights, through radial basis functions in this study. The efficacy of the proposed procedure is 
validated for a two-dimensional airfoil moving at different Mach numbers and pitching motions, under which 
the flow exhibits prominent nonlinear behaviors. All results confirm that our nonlinear model is accurate and 
stable for predicting not only aerodynamic forces but also detailed flowfields. Moreover, the model is 
robustness — accurate for inputs considerably different from the base trajectory in form and magnitude. This 
modeling preserves nonlinearity of the problems considered in a rather simple and accurate manner. 


I. Introduction 

A erodynamics problems are highly nonlinear and are routinely solved nowadays by computational fluid 
dynamics methods through validated computer codes, thanks to the availability of fast-improving computer 
technologies. However, the need for striving to reduce the computation time has become as strong, if not more, as 
ever because more details (thus more grid points and higher fidelity) are demanded and coupling of disciplines 
(such as structure, control, chemistry, etc.) required for considering an engineering system. Furthermore, 
optimization problems require repeated computations of similar setup, but with variations in design variables that 
may have a large variation in values. Thus, there exists a strong desire to replace the original full model with an 
approximate model that is much computational intensive. This may be achieved in several ways, such as model 
order reduction whose concept is to find an approximation with number of unknowns much smaller than the full 
system. Or one can base on a set of available data (empirical or CFD solutions), valid at discrete parametric 
conditions, to build a model that presumably will be valid inside the parameters (state) range but also beyond to 
some extent. The latter is typically employed in surrogate modeling, metamodeling, and the like. Notwithstanding, 
the key question (hence criterion) is that the fidelity of the original model must be preserved to the extent 
necessary. This requirement broadly translates into the preservation of nonlinear behavior in aerodynamics. This 
is the goal of the research reported here: modeling for efficiency and accuracy. 

Our research into model reduction has been motivated by the need for performing multidisciplinary analysis 
and optimization for an aeronautical system under the current Fundamental Aeronautics Program. In particular, 
aeroelasticity is the subject of interest — we are concerned with the coupling of unsteady aerodynamics and 
structure dynamics and its interactive effects, specially the limit cycle oscillations (LCO). For which, numerous 
scenarios must be evaluated to relate the output in response to the input parameters. 

Conventional aeroelastic stability analysis is often formulated simply in terms of an input-output relationship, 
where input is aerodynamic force and output is structural response. System identification (ID) method is a 
straightforward framework to construct a simple yet theoretically rigorous reduced order model (ROM) at the 
expense of detailed flowfield information. For example, such flow flied characteristic will be required to find a 
new aerodynamic shape to minimize occurrence of structure flutter in an aeroelastic optimization problem. 
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To retain flowfield details, the state-space project method based on modes superimposition has often been 
employed. Unfortunately, the linear superimposition of modes limits the method from describing nonlinear 
phenomenon such as limit cycle oscillations. 

The basic concept of State-Space Projection method is to find a subspace or transformation matrix, onto which 
the full order system can be projected with a much smaller dimension. A well-known method to define such 
subspace is proper orthogonal decomposition (POD). Beran et al. [2] have used the POD method to construct 
ROM of aeroelastic system. A key ingredient of POD is to collect snapshots of the time-dependent solution. Once 
the snapshot matrix is formed, the subspace can be derived from an eigenvalue analysis. Traditional POD method 
merely takes input characteristics into account. Moore [3] introduced the Balance Truncation (BT) method to 
define the subspace by including both the input and output characteristics of the full order system. As a result, the 
BT method often yields not only more accurate but also lower order ROM than POD. Laub [4] proposed a more 
efficient way to derive the transformation vector. Willcox and Peraire [5] and Rowely [6] respectively applied 
POD snapshot ideas to extend the original BT method to study fluid dynamics problems. 

Linear or nonlinear full order model representation is required for current POD/snapshot-BT method. It is 
relatively straightforward to obtain a linear model of a highly nonlinear system by retaining the first-order terms; 
this is simply done, for example, by finite differencing a CFD code. However, it remains challenging to obtain a 
nonlinear model of the CFD code, as elaborated below. Let us consider a time-dependent nonlinear system, 

x = f(x,u) ( 1 ) 

where x is a N*\ state vector; u a Pxl input vector. An appropriate projection can be chosen to construct a 
nonlinear ROM in a subspace <f > , expressed as: 

*,= 07 ( 0 * >«) (2) 

where x r is a R*\ state vector, with R « N. An intuitive way of arriving at a nonlinear representation is to 
include the second-order terms in the Taylor’s series expansion, thus involving the Hessian matrix of f However, 
it renders computationally intractable to obtain the Hessian matrix for a fluid system in CFD. 

Rewienski and Whitefl] introduced a trajectory piecewise-linear (TPWL) approach for constructing a 
nonlinear model: including the projections of state variables and the governing nonlinear differential equations 
onto a subspace and a proper weighting of a number of these truncated linear models expanded about some 
operating points (or called sample states in this paper as in [7]) along the state trajectory. Gu and Roychowdhury 
[7] proposed a similar but distinguishably different technique and demonstrated its advantages over the TPWL. 
Gratton and Wilcox [8] combined the TPWL and POD methods to derive a reduced model to study a flow control 
scheme for supersonic diffuser. Jesmani et al [9], taking advantage of linear characteristics of flow in streamline 
coordinates, applied the TPWL for prediction of pressure and water saturation in oil reservoirs. 

Previously, we presented the idea of assembling a number of weighted neurons associated with different 
inputs to arrive at a nonlinear reduced model via recurrent artificial neural network training [10] in which the 
weights are based on the radial basis functions. This thinking follows the general concept of metamodeling. The 
present work extends the similar concept for constructing a nonlinear model by collecting a set of submodels that 
are each valid in the neighborhood of sample points. The idea of TPWL comes handy so that we can simply 
replace the neurons with the piecewise linear solutions valid around the sample points. It is expected that 
including the first-order Taylor expansion terms will enhance the “ball” of validity, consequently the resulting 
model should be more robust. 

The objective of our work is to ultimately establish an accurate, stable, and efficient procedure to construct a 
nonlinear model that can be used reliably to describe nonlinear behavior in an aeroelasticity system; this paper 
presents the first installment towards that objective. Major efficiency gains will come in an upcoming paper in 
which the balanced truncation technique is implemented to reduce the order of the approximate system. The 
present paper is organized as follows. First, the basic concept of piecewise linearization (PWL) in a state 
trajectory is presented for a simple nonlinear scalar equation to illustrate how the PWL-based approximate model 
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works. Section III gives the details for constructing a nonlinear model from CFD solutions. Validation of the so- 
constructed model is shown in Section V and VI, using the AGARD test cases, CT5 and CT2 [13], as a baseline. 
Key issues to its construction are discussed in Section VI, with additional results attesting the efficacy of the 
proposed modeling approach. Finally, we summarize the paper and outline future plan for continuing the work. 

II. Piecewise-Linearization 


In this section the basic idea of piecewise linear approximation is outlined, while the interested reader should 
consult [1] for details. Applying the first order Taylor expansion of the nonlinear system Eq. (1) around sample 
points (x u Mj), (* 2 , u 2 ) ... (x p , u p ) gives 
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Each of the above expansion is assumed to have an adequate range of validity in the neighborhood of its 
sample point. It is anticipated that a proper assemblage of these piecewise linear approximation will provide an 
approximation over an enlarged domain and a nonlinear combination of them will potentially preserve the 


nonlinearity of the original nonlinear system. Introducing 
approximation near u,) to form a weighted sum of Eq. (3), 
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Eq. (5) can be rewritten, using the Jacobians of / with respect to state variable * and input variable «, as 
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( 6 ) 


This provides us a framework to represent the nonlinear model Eq. (1) with a combination of linear model, with 
three questions to be answered: (1) How many sample points (p) are necessary? (2) How are the sample points 
distributed? and (3) How are the weight functions determined? We are not aware of any rigorous and general rules 
addressing these issues and suspect that one may likely rely on experiences on practicing and verification and need 
to have good understanding of the physical problems under consideration. Nevertheless, some rules and 
knowledge, albeit problem dependent, will be accumulated, thus providing a useful guide. The first two questions 
will be evaluated quantitatively in the present study. The third one has been studied in the general context of 
surrogate model [11]. Here, we shall adopt the radial basis function (RBF) of the Euclid distance between the state 
* and the sample (operating) point, also used in [10], 
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where || x-x t || is the L2-norm (Euclidean distance) of the vector jc-jc„ jc, is the sample state, Si is the width parameter. 
In this paper we take S k =1/V50 for all CFD calculations. Schematically the nonlinear model can be viewed as in 
Figure 1. 

linear module 



Figure 1. Assembly of piecewise linear module to form an approximate nonlinear system. General Guidelines. 


A. Validation of the Piecewise Linear (PWL) Model 


In order to evaluate the validity of the above piecewise linear concept for a nonlinear system, a simple model is 
considered, 

x = /(x,w) = 5-sinx + x + w ^ 

According to Eq. (4), the nonlinear system Eq. (8) is approximated as, 

p 

x - /(x,w) = 5*sinx + x + u ~ w. *(/ 0 + tf. (x-x.) + Z>. [u- 

w ' 1 * (9) 

Here, the notation a, and bi is used to replace Jacobian A h B iy since the number of states of the model is one. These 
derivatives a, and bi are available analytically 
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( 10 ) 

Obviously, accuracy of Eq. (9) highly depends on the choice of sample points ( jc „ u,) at which the training of 
input-output relationship will be imbedded in the model. The training signal can be determined according to the 
range of interest. Here, the input signal is defined as, 

u = sin(6r) + cos(6/) 

With the initial condition x\ = 0 and input signal Eq. (11), the output of Eq. (8) can be numerically integrated 
by using the standard 4-stage Runge-Kutta method for time integration and is shown in Figure 2. 


A total of 30 sample states x were chosen uniformly to build the piecewise linear model, Eq. (9), with the first 
and last sample points set at t = 0 , and t =3.0 respectively. Figure 3. is the output of the PWL model with different 
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numbers of PWL models: one PWL is only valid around the first sample point as expected. Using the error 
between the “exact” output ( x ) and the PWL model 


Error = 


1 N 

lOOx — V \x. -x.\ 
' f| 


allows us to evaluate the effects of the size {p ) of PWL models. 


( 12 ) 


As shown in Table. 1, the error decreases as the number of PWL increases. Eventually, a 3-PWL model was 
chosen to form a representation of the nonlinear system Eq. (8) since p = 5 offers only an additional 0.12% 
reduction in error. As stated above, the PWL was obtained by training the nonlinear system Eq. (8) with the signal 
Eq. (11) for tE[0, 3.0]. To see the usefulness (or robustness) of this approximation, it is necessary to evaluate the 
model’s performance at t > 3.0(s) or with different signal input. Figure 4. is the output using the 3-PWL derived 
before, in comparison with the solution of the full nonlinear system. The 3-PWL remains accurate for t > 3.0(s). 
This is not entirely unexpected because the state variation after t > 3.0(s) still falls in the range of x(t) for t< 3.0, 
but nonetheless confirming that the model also adjusts to the asymptotic periodic solution after the transient period. 
Next, we test a very different input signal 

u = sin’(6/) (13) 

Figure 5. shows that the same 3-PWL still gives an accurate output. However, it is worth reminding that the 
reliability of the model strongly depends on the training data, namely the sample points; some judicious choice of 
training trajectory is useful for maintaining the model’s range of validity. 



time(s) 

Figure 2. Output of system Eq. (8). 



time(s) 

Figure 3. Output of PWL model. 


Table. 1 Error of different PWL No. 


PWL No. 

Error (%) 

2 

8.69 

3 

1.19 

5 

1.07 
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Figure 4. output at t > 3.0(s). 




(a) Input signal. 


time(s) 

(b) Output comparison. 


Figure 5. results with input signal (13). 


III. CFD modeling by Piecewise Linear Construction 


A 2D inviscid CFD solver is used here to demonstrate how to construct a CFD-based PWL model. From the 
aeroelasticity viewpoint the output of CFD is aerodynamic forces, which often exhibit nonlinear behavior with 
respect to inputs, such as motion or shape change of a wing. Without loss of generality, the inviscid CFD equations 
can be written as. 


a(uq) 

dt 


R(q,u,v) 


(14) 


Where V is the volume of a computational cell, q consists of conservative flow variables, R(q,u,v) is the residual 
describing derivatives of flow flux functions, for a set of input conditions, u and v, which are the structure 
displacement and velocity respectively for the cases considered in this paper. The spatial discretization contained in 
the residual term R(Q) can be determined by a state-of-the-art numerical flux method; here we employ the AUSM + - 
up scheme [12] for the flux function together with the van Albada limiter in the MUSCL interpolation of primitive 
variables. The time integration is done by using the 4-stage Runge-Kutta method together with a dual-time iterative 
scheme incorporated to ensure that the accuracy of the time-discrete form of unsteady equations be preserved at each 
time step. 


Now for the model system, we follow the same procedure described in Eqs. (3-5). The residual term R(q,u,v) in 
Eq. (14) may be approximated with a collection of the first-order Taylor series terms expanded about a set of 
appropriately defined sample points, Qr= (q„ U/, v,-), i =1,..., p. Expanding R(q,u,v) about a sample state Q„ we get 


R(Q) = R(Q,)+A 1 (q-q,) + B u («-«,)+B !J (v-y,) + ... 
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where the Jacobians of R(q,u,v) with respect to its arguments are A, Bj and B 2 respectively. Hence, Eq. (14) 
becomes, 
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(16) 
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The right-hand term of Eq. (16) is a weighted combination of locally linearized models; it is completely 
linearized except the weighting functions. The solution of this system is computationally tractable, bypassing the 
need of flux evaluations for the entire flow domain and requiring only matrix-vector multiplications besides 
evaluating the weights. In this study, the same procedure used in the full model for time integration is again used 
here. It is worth noting that in this formulation, the grid motion expressed in terms of v need be computed as well, 
via either a prescribed motion or an implicitly coupled system as in the case of flexible wing whose deformation is 
affected by and changes the flow simultaneously. 


Assuming a continuous grid motion, the time derivative can be further rearranged, by knowing the volume v is a 
function of airfoil displacement u and velocity v = du/dt , 




dq 


3q 


dt' 


du 


dt 


dt 


Substitute Eq. (17) into Eq. (16), 

I7 = X^( R (Qi)+ A , (q -q,)+ B u ( u - ' u /) ■ + B 2,/ ( v - v ,) - c, v) 


(17) 


(18) 


Now the grid motion is explicitly folded in the right-hand-side term through the Jacobian matrix C; all that is 
required boils down to the time evolution of the state vector q. As we consider only rigid body motion in this study, 
this is the equation used. Thus, the following dual-time iterative system is solved for the increment of q in each time 
increment of the 4-stage Runge-Kutta method, 


S e ^ + S a ^,^i( R(Q()+A| ( q ., | ) +Bu („_ u J + B !j ( v — Vi )-c,v) 

/-I 1 


(19) 


where the subscript “m” denotes the iterative index for the dual time and subscript “n” for the physical time. 

The fictitious step used in the dual time At is calculated by, 

Ar = ]?w l At loco/l 

i=l ( 20 ) 
Solving Eq. (19) does not require tracking of grid deformation, the grid motion is already included in Qv. 
However, aerodynamic force is still obtained by integrating the pressure force over the instantaneous position of the 
moving body. The grid information should be updated at each physical time step. For a 2D problem, like the cases 
studied in this paper, the computational time for grid deformation is negligible and so is the integration of forces 
over the body. However, by skipping a direct integration of forces over the body at each time step, it is consistent, to 
the order of accuracy involved in the above weighted combination of piecewise linear expansions, to express the 
aerodynamic forces in the form of the first-order Taylor expansion of the aerodynamic force about each sample 
point, 


F-F. + |S 
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Then, the PWL model for the aerodynamic force is formulated similarly using the same set of sample points 
and weighting functions. 
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(22) 



All three procedures have been implemented to compute the aerodynamic forces and they do not give rise to 
noticeable differences, hence not shown in the present paper. The aerodynamic force shown in the remainder of 
the paper is computed through Eq. (22). 


IV. Validation of Nonlinear CFD Model 

In order to demonstrate the validity and performance of the above-constructed CFD PWL-based nonlinear 
model, a 2D NACA 0012 airfoil in pitching motion will be used for various flow conditions that result in vastly 
different aerodynamic behaviors. Special interest is paid to the nonlinear characteristics. In the present study, we are 
exclusively interested in unsteady aerodynamic characteristics where the inviscid mechanism plays the dominant 
role. In order to confirm the baseline solution, some standard test cases put forth in the AGARD report [13] will be 
used, specifically the so-called CT2 and CT5 cases. These two cases have different flight Mach numbers and 
pitching parameters, as summarized in Table 2. The pitching motion of the airfoil is described by its angle of attack 
changing in time 


a{t) = a m + a Q sin (cot) = a m + a 0 sin(2 Id), 




(23) 


where k is the reduced frequency defined by the free stream velocity U*, and the chord length c, and the pitching 
centers at x m , given in Table.2 as well. 


Table. 2 Parameters of CT2 and CT5 cases 


Case 



a 0 

k 

Xm 

CT2 

0.6 

3.16° 

4.59 

0.0811 

0.273 

CT5 

0.755 

0.016° 

2.51 

0.0814 

0.25 


Since the inviscid force is the primary controlling mechanism for the unsteady aerodynamic characteristics of 
interest, an O-type grid around the airfoil is appropriate for the CFD computation. A grid of 61x21 points (in radial 
and circumferential directions respectively) is generated, as shown in 

Figure 6. This results in a state vector q with a size of 4800 (60x20x4, where a CFD solution vector for a 2D 
problem has 4 components [p, pu, pv, pE] for each computational cell.) 
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First, we validated the accuracy of the CFD procedure against the experimental data reported in [13]; the 
comparison of aerodynamic forces in a cycle for both the CT2 and CT5 cases is given in Figure 7. and Figure 8. 
respectively, showing the CFD results in good agreement with the data. More importantly, both test cases exhibit a 
strong nonlinear behavior in the moment coefficient and nonsymmetrical forces between downward and upward 
motions. It is noted that CT2 has a significantly larger excursion in pitching than CT5, thus it is expected to show 
stronger nonlinear effects. On the other hand, CT5 operates at a higher Mach number, and it is anticipated to have a 
supersonic pocket in the flowfield during certain phase of the pitching. 

This validation establishes a solid foundation for trusting the solutions at any point along the pitching trajectory; 
that is, any point is equally legitimate to be chosen, insofar as its accuracy is concerned, as a sample solution-a 
piecewise linear model in this case. In the remainder of this paper, we shall focus on the performance of the 
nonlinear model built from such PWL solutions (submodels). 




a a 

Figure 7. Validation of CFD solution for the AGARD CT5 problem. 




Figure 8. Validation of CFD solution for the AGARD CT2 problem. 


V. Practical Issues about Constructing the Nonlinear Model 


As mentioned previously, the training dataset, i.e. the piecewise linear solutions, must be representative, 
containing essential flow characteristics, so that the final overall model can be useful for a range sufficiently far 
away from the baseline trajectory to cover a wide range of input parameters. The number of sample points should 
not be too “large” so that the computational effort is practical and the distribution of these points should not be too 
“specific” so that no a priori knowledge about the flow details is necessary. However, some intuitive judicious 
judgment may be helpful. 

A. Effect of Sampling 
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To our knowledge, there are no intelligent and general criteria to pre-determine a better set of sample points. We 
simply base on the time coordinate, instead of following the flight trajectory of angle of attack. During a pitching 
period T, the physical time step for time integration is divided into 60 equal time steps, i.e., At = T/60. Out of these 
60 time slices, 20 equally-spaced sample points are chosen, at which the PWL solutions, hence denoted as PWL(20), 
are gathered, as shown in Figure 9. for the CT5 problem. For constructing the CT5 model, we used a 0 = 3.5° as the 
baseline training. To test the effect of the number of PWL solution, we also used 13 sample points, clustered around 
the extreme points of the motion; the solutions, called PWL(13), are included in Figure 10. for comparison with the 
full CFD and the PWL20 solutions; again both the PWL-based solutions agree closely with the CFD solution, as 
shown in 

Figure 11, also for other two flight trajectories with a 0 = 3.0° and 2.5°. There does not appear a definitive trend to 
judge whether PWL(20) performs better than PWL(13), however the latter certainly requires much less 
computational efforts. As a result, this does suggest a more intelligent way (in this case clustering sample points 
near the peak angles) can be more effective. 



Figure 9. PWL(20) sample points. 


input 



Figure 10. PWL(13) sample points. 




a 

(a) a 0 =3.5°. 


a 
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a a 

(b) a 0 =3.0°. 




a a 

(c) a 0 =2.5°. 

Figure 11. Effect of distribution of sample points on CT5 solutions: nonuniformly-spaced (PWL(13)) vs. 
uniformly spaced (PWL(20» points. 

Shown in 

Figure 12. are the results of the full CFD model and the nonlinear model using 30, 20 and 10 PWL solutions for 
CT5 at a 0 = 3.5°, suggesting that both 20 and 10 PWL models can faithfully capture nonlinear features of the 
original CFD system as well. On the other hand, the linear results, obtained by retaining only the l sl order Taylor 
expansion around the steady operation point, which is at the time t = 0(s), completely misses the nonlinear 
phenomenon in the original system. It is noted that the linear model is commonly used to generate linear model of 
CFD system for order reduction. As the linear model does not contain any time varying information in the model, it 
becomes clear collecting time-varying submodels into the model is a very essential element and a key to maintaining 
accuracy of the model. After the training, the so-built PWL-based model can be used for other conditions. 

Figure 12 (b) and (c) are the results at a 0 = 3.0° and 2.5° respectively. Not surprisingly, the PWL-model is also 
valid, since both conditions fall within the training data range. The linear model results are also included for 
comparison, again unable to capture nonlinear phenomenon at lower angle of attacks. 




11 

American Institute of Aeronautics and Astronautics 


(a) a 0 =3.5°. 

Euler 




a a 

(b)a„=3.0°. 

Euler 




a a 

(c) a 0 =2.5°. 

Figure 12. Effects on accuracy by the number of PWL solutions used. 

Next there arises a question as to how the model trained for the condition at ao= 3.5° and k = 0.0814 behaves at 
different conditions. First, we change the reduced frequency to k- 0.1, a 20% larger than the baseline value. Shown 
in Figure 13, the PWL-model still holds valid. The 10-PWL solution for this case begins to show some noticeable 
departures from the full solution at the peak Cm, indicating more sample points are needed. It is noted that even 
though the pitching frequency does not explicitly appear in the formulation, it implicitly affects the flow solution q, 
which in turn changes the residual R and related quantities. Obviously, if pitching angle or frequency changes flow 
structure significantly from those at sample points, then the accuracy of the PWL-model will deteriorate. 




a a 

Figure 13. PWL results at different pitching frequency k = 0.1 (a 0 =3.5°). 

It is also necessary to check if the construction of the solution procedure guarantees convergence. This is tested 
by decreasing the time step size, in this case by halving the physical time step, dt = 77120; results in Figure 14. 
confirms convergence. 
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a a 

Figure 14. PWL results using a half time step dt = T/120 (a 0 =3.5°). 

Now we consider the scenario of having an input, the flight motion, considerably different from the baseline. In 
this case we impose a cosine signal with the other parameters unchanged. 

Figure 15. shows that the PWL-model is still valid with a different input. 




a a 

Figure 15. PWL results with cosine signal (a 0 =3.5°, dt = T/60). 

With the integrated quantities thoroughly compared, it validates the accuracy of the PWL-based model. It is also 
interesting to compare the detailed flow structure against the full CFD solution, as shown in Figure 16, and again it 
confirms the reliability of the present nonlinear model. 



Euler PWL 

t=T/4 
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Euler 


t=T/2 


PWL 


Figure 16. Flow-field expressed by pressure contours at two time slices, cosine signal (a 0 =3.5°, dt = T/60). 

Finally, we compare the computational cost for performing a full order CFD analysis and the additional work 
required to construct the PWL-based nonlinear model, as shown in Figure 17. One should note that this extra 
overhead is one time only; any computation for other conditions will require a very minimal effort to perform 
matric-vector multiplications and time integrations, while the full CFD model will be another brand new 
operations. Hence, in this sense, the so-constructed nonlinear model is not only accurate but also useful for 
parametric studies and design optimization. 
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PWL Euler calculation 

Figure 17. Time consumptions in the full Euler solution and PWL-based modeling (Intel core CPU i7). 

B. Further Application of the PWL-based Modeling to CT2 Problem 

In this section, we further apply the previously established PWL-based nonlinear model for another unsteady 
pitching airfoil problem, CT2 in [13]; this problem, while remaining completely subsonic during the entire cycle, is 
more difficult aerodynamically than CT5 in the sense that the aerodynamic forces display a stronger nonlinearity, as 
seen in Fig. 8, where the lift force has a noticeable nonlinear behavior near high angle of attack, exhibiting 
asymmetry in high angles of attack. It is interesting that even at this high angle, the inviscid mechanism is still 
dominant as validated against the data. The training is carried out at a 0 = 8°, again using 20 PWL solutions. 

In Figure 18, we show the performance of the PWL-based solutions against the full Euler solutions of different 
pitching amplitudes. At a small amplitude (a 0 = 1°), the flow behaves linearly and the linear model is close to the 
Euler solution. On the other hand, when the amplitude is only moderate, nonlinearity first appears in the moment 
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coefficient; then the lift coefficient exhibits a strong nonlinear phenomenon with a double loop, suggesting a second 
higher-frequency flow structure is embedded in the main flow. The linear model clearly is incapable of following 
the correct flow behavior. 







a a 

(c) a 0 = 8°. 

Figure 18. Comparison of Euler, nonlinear PWL-based model and linear solutions of the CT2 airfoil at 
different pitching amplitudes. 

A snapshot of the detailed flowfield is displayed in Figure 19. for a high angle motion (a 0 = 8°), showing no 
discemable differences between the Euler and PWL-based solutions. 
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Euler PWL 


t=T/4 



Euler PWL 

t=T/2 

Figure 19. Pressure contours of AGARD CT2 at two time instants (T/4 and T/2) for the airfoil pitching at a 0 

= 8 °. 

Now it is extremely informative to see how the model we constructed performs outside the training trajectory. 
With the previously-trained model for the trajectory prescribed by Eq. (23) with ao=8°, and a m =3.16°, we applied it 
to a new motion defined by, 

a= a m + a 0 sin 3 (2 kt) (24) 

with a 0 = 8°, 0^=3. 16°, its flight trajectory shown in Figure 20, contains a high (triple) frequency content resulting 
from the cubic function. 
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Figure 20. New input trajectory defined in Eq. (24). 

A vastly different behavior is revealed in both the lift and moment coefficients, Figure 21 showing sharp 
discontinuities at the mean angle a m = 3.16° in both upward and downward phases. The knot point (a ~ 6°) in Q 
does not find a particularly unusual behavior in C m . At the mean location, the pitching velocity nearly vanishes, 
unlike the CT2 motion Eq. (23). 




a a 

Figure 21. Lift and Moment coefficients responding to Eq. (24) input. 

Taking a notch higher by setting do = 10°, upward away from the training trajectory, we show the results 
predicted by the model in Figure 22. Comparing them against the full Euler solution, this case further confirms that 
the current PWL-based nonlinear model is capable of accurately predicting flowfields considerably beyond the 
training trajectory., except with some discrepancy near the peak angle. 




Figure 22. Lift and Moment coefficients responding to Eq. (24) input. 
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Finally we investigate the flowfield at the knot points, or at the time slices t=T, 3T/2, as shown in Figure 23; the 
results computed by the model are nearly identical to the CFD results in every detail of the flowfield, also indicating 
a nonsymmetry here even though the airfoil is at the same angle (a = 3.16°) at both times. 



t=T 



t=3T/2 


Figure 23. Detailed flowfield by the PWL-based model and CFD, a 0 = 10°. 


VI. Conclusion 

A nonlinear modeling method is proposed by nonlinearly assembling a number of local piecewise linear 
solutions obtained at some sample points. These points may be defined along the physical or state-variable 
trajectory, for example simply along the time coordinate in the problems considered in this paper. We show how the 
piecewise linearization is carried out for the CFD applications. The validity and efficacy of the PWL-based 
nonlinear model is confirmed first for a scalar equation and subsequently for AGARD test cases, CT2 and CT5; the 
results shows the accuracy of the model not only for the baseline conditions, but also for conditions beyond the 
training trajectory. The nonlinear characteristics in the flowfield and aerodynamic forces with respect to the varying 
angle of attack are predicted in close agreement with the full Euler model. The model is capable of predicting the 
flowfield accurately in response to an input that is vastly differently from the baseline, yielding a flow structure not 
existing in the baseline solution. The accuracy, robustness, and efficiency of the model constructed by the proposed 
method suggest that the method can be valuable for the situations where repeated computations are required for a 
host of input and parameters changes, such as design optimization. 
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To further reduce the computation cost for preparing the PWL-based model, some model reduction techniques 
by projection onto a much smaller subspace, such as balanced truncation, has been implemented and the results will 
be presented elsewhere. 
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